Para-hydrodynamics from weak surface scattering in ultraclean thin flakes

Electron hydrodynamics typically emerges in electron fluids with a high electron–electron collision rate. However, new experiments with thin flakes of WTe2 have revealed that other momentum-conserving scattering processes can replace the role of the electron–electron interaction, thereby leading to a novel, so-called para-hydrodynamic regime. Here, we develop the kinetic theory for para-hydrodynamic transport. To this end, we consider a ballistic electron gas in a thin three-dimensional sheet where the momentum-relaxing (lmr) and momentum-conserving (lmc) mean free paths are decreased due to boundary scattering from a rough surface. The resulting effective mean free path of the in-plane components of the electronic flow is then expressed in terms of microscopic parameters of the sheet boundaries, predicting that a para-hydrodynamic regime with lmr ≫ lmc emerges generically in ultraclean three-dimensional materials. Using our approach, we recover the transport properties of WTe2 in the para-hydrodynamic regime in good agreement with existing experiments.

In this supplementary information we present the calculation of the electron-electron mean free path, derive the operator of angular diffusion, and detail the solution of the differential equation for boundary scattering.

ELECTRON-ELECTRON COLLISION RATE FROM THE SCREENED COULOMB INTERACTION
To quantify the strength of the interaction we focus the study on the relaxation length, l ee , the mean free path for electron-electron collisions. In Fermi liquid theory normal metals have a relaxation length with a T −2 temperature dependence. We emphasize that in WTe 2 , the assumption of a T −2 -dependence is not justified a priori. Instead, it is necessary to calculate l ee down to low temperatures in order to extract a reliable temperature dependence. In our previous calculation [1], we indeed found that the low temperature limit approaches a T −2 behavior, but such happens only at energies much smaller than the Fermi level, which in WTe 2 is unusually low at E F ≈ 150K. Since the Fermi energy, and the Fermiology of the Fermi surface play a prominent role in WTe 2 [2], and since our previous estimate for ℓ ee ≈ 10µm was only 1.5 orders of magnitude to large to explain the hydrodynamic electron flow, here we repeat our calculation for two other band structure candidates for WTe 2 . This way, we can make sure that the reported absence of a sufficiently small ℓ ee is not an artifact of the specific choice of band structure, but holds invariably for this material.

Formalism
Electron-electron interactions arise from the coupling of the electrons to the electromagnetic field. In 3D, the bare coulomb interaction is given by U bare (q) = 4πe 2 |q| 2 . The presence of electrons modifies the bare interaction, such that the renormalized interaction W becomes . (S1) Here, Π(q, ω) is the particle-hole susceptibility of the electromagnetic field, given by the one-loop integral i, j are band indices, n(ϵ) is the Fermionic occupation number, ξ = ϵ − ϵ F and |ρ i,k,j,k+q | 2 = |⟨u i,k e iqr |u j,k+q ⟩| 2 is the exchange density. At finite temperature the electronic self-energy due to the Fock diagram is where b(ϵ) is the bosonic occupation number. The electron-electron relaxation rate of a state is defined in terms of the self energy as The average relaxation rate is then the average of the relaxation rate on the Fermi surface, normalized by the total density of states: Finally, the relaxation length is given as the inverse of the relaxation rate,

Candidate band structures
Weyl semimetal candidate band structure The chemical potential resides at µ = 7.22eV , the lattice parameters are a x = 0.3483nm, a y = 0.6265nm, a z = 1.4043nm. The resulting parameters are listed in Supplementary Table I. Relaxed candidate band structure The chemical potential sits at µ = 8.5083. This band structure is motivated by Ref. [2], with lattice parameters a x = 0.3460nm, a y = 0.6200nm, a z = 1.3090nm. Compared to the other band structures, this corresponds to a compression along the k z -direction. The position of the atoms inside the lattice was subsequently optimized to minimize the total energy. The result is shown in Supplementary Table II. Candidate band structure based on hybrid functionals Here, the chemical potential is µ = 6.965eV . The chosen lattice parameters are a x = 0.3483nm, a y = 0.6265nm, a z = 1.4043nm. Supplementary Table III lists the band structure parameters obtained when using hybrid functionals in the self consistent stage of the calculation in order to better account for exchange. We briefly comment on consistency checks using the density of states, and why they are less relevant in the present calculation. For example, the total density of states of the relaxed phase is 7.81eV −1 nm −3 and according to the calculation of the real part of the susceptibility at ω = 0.01eV it is 6.47eV −1 nm −3 , which agrees within 20% with the total density of states. Likewise, the total density of states of the Weyl semimetal phase is 7.12eV −1 nm −3 and according to the calculation of the real part of the susceptibility at ω = 0.01eV it is 4.27eV −1 nm −3 , which is further apart, with a 67% difference. As the main reason we suspect the more sensitive dependence on ω in the Weyl semimetal band structure. However, for a faithful estimate of the self-energy this effect is not relevant because we are predominantly interested in the imaginary part of U , which becomes zero exactly at ω = 0.

Results for the electron-electron mean free path
The electron-electron mean free path according to Supplementary Eq. (S6) and its temperature scaling are shown in Supplementary Fig. 1 for all three candidate band structures. Clearly, the low-energy asymptotic scaling onsets at or even slightly above 150K in all three cases. Extrapolating the obtained temperature scaling, we arrive at the low-temperature quoted in the main text.
Supplementary Figure 1. lee as function of temperature for the three band structures. (a) Small pockets: due to the compensation of electrons and holes and the similar velocities of the bands, lee and l hh are similar, we can see a good fit to a power law, but with a weaker dependence on temperature than inverse square. (b) Anisotropic pockets: the Fermi level was chosen in a way that the compensation of holes and electrons is not perfect, more electrons than holes. This is reflected both in an effective lower Fermi level and lower velocity for the holes, which manifest in the form of l hh < lee. Still, we can see a good fit to a power law, specifically for electrons, a T −2 . (b) Relaxed phase: resembels a metal and as such fits well to a T −2 power law. The relatively large lee compared to the previous two band structures is mainly due to the large Fermi level, which leads to high screening, reducing the effective interaction.

DIFFERENTIAL OPERATOR
The derivation of the differential operator was first done in Ref. [3], but uses a highly customized notation. For this reason, we reiterate here the steps needed to obtainÔ θ . We start from the general boundary scattering condition in terms of the scattering potential (Eq. (2) in the main text), f > (k) = f < (k) + I where the integral is defined as Recalling that the scattering potential has the explicit form W (k ′ , k) = W 0 e − |k−k ′ | 2 2σ 2 , we proceed with a saddle point approximation. Denoting g(k ′ ) = k ′ z (f < (k ′ ) − f < (k)), it is clear that g(k) = 0. Expanding up to second order in k ′ − k, in spherical coordinates on the Fermi surface, it is therefore Inserting the expansion into the integral I = I 1 + I 2 , we obtain In polar coordinates relative to k, we notice that (k ′ − k) i is anti-symmetric with respect to the transformation is symmetric. Therefore the integral on the components that are perpendicular to k, meaning i = θ, φ vanish. On the other hand, the integral on the radial direction, does not vanish from this argument, but instead vanishes identically, as |k ′ | = |k| is enforced identically on the entire integration range. Thus it is I 1 = 0.
For I 2 , obviously the same argument holds for the radial direction. It thus suffices to look at directions i, j that are perpendicular to k. Considering i ̸ = j next, the integrand is odd in either direction and vanishes due to antisymmetric cancellation. So the only terms that contribute are i = j for the two directions perpendicular to k, The integral inside the brackets can be solved by noticing that for the two directions perpendicular to k, the result is the same from rotational symmetry around the k/|k| axis: In summary, the integral I evaluates to We note that ∂ 2 i where i = θ, ϕ is exactly the Laplace-Beltrami operator on the sphere. We can therefore write more formally (S14) Using the identity ∇ 2 (k z f < (k)) = f < (k)∇ 2 k z + k z ∇ 2 f < (k) + 2∇k z · ∇f < (k) and noting that ∇ 2 k z = 0, we recover the form of Ref. [3], In spherical coordinates, this is Since we consider a thin slab where the thickness d is much smaller than the width w, we can neglect derivatives with respect to the in-plane angle ϕ, which are expected to only contribute weakly due to the smooth dependence of f on ϕ.
In term of the reduced distribution function h(z, θ) the boundary condition therefore simplifies to h(d/2, −|θ|) = 1 + Q cos θ cot θ∂ θ (sin θ∂ θ ) − 2 sin θ∂ θ h(d/2, |θ|) (S18) which, after shifting θ → θ − π/2 and adding the specularity coefficient R θ is the form quoted in the main text in Eq. (3). While it may look as though Supplementary Eq. (S19) contains a divergence at at θ = 0, we point out that this divergence is of the form 1/ sin θ, which is rendered finite upon integration with the Jacobian k 2 F sin θ on the sphere.

SOLUTION OF THE DIFFERENTIAL EQUATION FOR BOUNDARY SCATTERING
The differential boundary condition, Eq. (5) in the main text, can be rewritten in terms of the variable s = sin θ, and reads explicitly A unique solution for c(θ) is obtained by imposing smooth continuity for this periodic and even function, i. e. we demand that c ′ (0) = c ′ (π/2) = 0. However, since θ = 0, π/2 are singular points, this way of solving Supplementary Eq. (S20) is impractical, and it is preferable to instead use initial conditions for c(θ 0 ) and c ′ (θ 0 ), where θ 0 ≈ 0 is chosen small enough for the effect of the derivatives in Supplementary Eq. (S20) to become vanishingly small at θ 0 . This makes sense because for θ → 0 Eq. (5) of the main text can be solved algebraically, yielding Then, using initial values c(θ 0 ) = 1 − (1 + λ)(1 − c 0 (θ)) and c ′ (θ 0 ) = (1 + λ)c ′ 0 (θ), and slowly adjusting the tuning parameter λ ≈ 0 such that the solution does not diverge at π/2 presents a reliable way to construct physically relevant solutions for c(θ). A typical iterative solution is presented in Supplementary Fig. 2 Unfortunately, this procedure still encounters issues at small α as the distribution function becomes exponentially close to 1, which means that the initial condition c(θ 0 ) ≈ 1 is rendered increasingly sensitive to machine-precision limitations. We found two avenues to remedy this issue.
One approach amounts to choosing a much larger θ 0 , which obviously cannot capture the form of c(θ) for angles below θ 0 , but converges smoothly to the correct solution at large θ. The second option is to expand Supplementary Eq. (S20) in the limit α → 0. Doing so removes the exponential functions, and thus the essential singularities at θ = 0, so that the choice of very small θ 0 remains unproblematic even for small α. Supplementary Fig. 2 shows the differences and commonalities between both approximations.